'''
WARNING:
  Preserve Python3 compatibility
  print()
'''


'''
Jorg Stelling model from
A synthetic tunable mammalian oscillator
'''

import sys
from numpy import *
from scipy import *
from bricks import *


class TristableSwitch:
  '''
  ARMUMENTS
    TIME : simulation time (minutes)
  '''
  gDEBUG = True
  gTIME = 10000
  gG1 = 200.0
  gG2 = 200.0
  gG3 = 200.0
  gG4 = 200.0
  gG5 = 200.0
  gG6 = 200.0
  gG7 = 200.0
  gpc_start = [];
  gpc_half_life = 0;
  grm_start = [];
  grm_half_life = 0;

  '''
  MODEL PARAMETERS

  G1: [57][pmin][TAL97:KRAB][TAL57:VP16][NicTAL:KRAB] [NEPTUN]
  G2: [97][pmin][TAL57:KRAB][TAL97:VP16][NicTAL:KRAB] [mCITRIN]
  G3: [Nic][pmin][TAL57:KRAB][TAL97:KRAB][NicTAL:VP16][RFP]
  G4: [12xZF-BS][P(PIR)ON][pmin][TAL57:KRAB][TAL97:KRAB][NicTAL:VP16][TAL95:KRAB]
  G5: [95][pSV40][p(PIR)ON][TAL97:KRAB][TAL57:VP16][NicTAL:KRAB]
  G6: [95][12xZF-BS][pIL-2][TAL57:KRAB][TAL97:VP16][NicTAL:KRAB]
  G7: [95][pCMV][PIP-KRAB][HetAct]
  '''

  


  def __init__(self, dbg, t, g1, g2, g3,g4,g5,g6,g7, pcs, pchl, rms, rmhl):
    self.gDEBUG = dbg
    self.gTIME = t
    self.gG1 = g1
    self.gG2 = g2
    self.gG3 = g3
    self.gG4 = g4
    self.gG5 = g5
    self.gG6 = g6
    self.gG7 = g7
    self.gpc_start = pcs;
    self.gpc_half_life = pchl;
    self.grm_start = rms;
    self.grm_half_life = rmhl;



  #################
  # SIMULATION PART
  #################



  def plotResults(self,TIME, T97KRABm, T97KRABp, T57KRABm, T57KRABp, T97VP16m, T97VP16p, T57VP16m, T57VP16p,  TNICKRABm, TNICKRABp, TNICVP16m, TNICVP16p, T95KRABm, T95KRABp, PBSm, PBSp, PBSa, HACTm, HACTp, HACTa, PC, RM, NEPm, NEPp, MCTm, MCTp, RFPm, RFPp):
    if sys.version_info[0] > 2:
      return

    from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
    from matplotlib.figure import Figure
    from matplotlib.dates import DateFormatter
    from random import random
    
    NONCESIZE = 1000;
    NONCE = int(random()*NONCESIZE)

    if self.gDEBUG:
      f0 = figure()
      f1 = figure()
      f2 = figure()
      f3 = figure()
      f4 = figure()
      
    else:
      f0 = Figure()
      f1 = Figure()
      f2 = Figure()
      f3 = Figure()
      f4 = Figure()
      
    f0.set_facecolor('white')
    s0=f0.add_subplot(111)
    s0.plot(TIME, HACTm, linestyle='dashed', label='RIP mRNA', color='blue')
    s0.plot(TIME, HACTp, label='RIP', color='blue')
    s0.plot(TIME, PBSm, linestyle='dashed', label='PIP mRNA', color='green')
    s0.plot(TIME, PBSp, label='PIP', color='green')
    s0.set_xlabel('Time [min]')
    s0.set_ylabel('Concentration [nM]')
    s0.legend(loc="upper right")
    lg = s0.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)

    f1.set_facecolor('white')
    s1=f1.add_subplot(111)
    s1.plot(TIME, PBSa, label='Pristinamycin-PIP', color='blue')
    s1.plot(TIME, HACTa, label='Rapamycin-RIP', color='red')
    s1.plot(TIME, PC, label='Pristinamycin', color='blue', linestyle='dotted')  
    s1.plot(TIME, RM, label='Rapamycin', color='red', linestyle='dotted')
    s1.set_xlabel('Time [min]')
    s1.set_ylabel('Concentration [nM]')
    s1.legend(loc="upper right")
    lg = s1.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
     
    f2.set_facecolor('white')
    s2=f2.add_subplot(111)
    s2.plot(TIME, T97KRABm, linestyle='dashed', label='TAL97:KRAB mRNA', color='blue')
    s2.plot(TIME, T97KRABp, label='TAL97:KRAB', color='blue')
    s2.plot(TIME, T57KRABm, linestyle='dashed', label='TAL57:KRAB mRNA', color='green')
    s2.plot(TIME, T57KRABp, label='TAL57:KRAB', color='green')
    s2.plot(TIME, TNICKRABm, linestyle='dashed', label='NICTAL:KRAB mRNA', color='red')
    s2.plot(TIME, TNICKRABp, label='NICTAL:KRAB', color='red')
    s2.plot(TIME, T95KRABm, linestyle='dashed', label='TAL95:KRAB mRNA', color='pink')
    s2.plot(TIME, T95KRABp, label='TAL95:KRAB', color='pink')
    s2.set_xlabel('Time [min]')
    s2.set_ylabel('Concentration [nM]')
    s2.legend(loc="upper right")
    lg = s2.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    
    f3.set_facecolor('white')
    s3=f3.add_subplot(111)
    s3.plot(TIME, T97VP16m, linestyle='dashed', label='TAL97:VP16 mRNA', color='blue')
    s3.plot(TIME, T97VP16p, label='TAL97:VP16', color='blue')
    s3.plot(TIME, T57VP16m, linestyle='dashed', label='TAL57:VP16 mRNA', color='green')
    s3.plot(TIME, T57VP16p, label='TAL57', color='green')
    s3.plot(TIME, TNICVP16m, linestyle='dashed', label='NICTAL:VP16 mRNA', color='red')
    s3.plot(TIME, TNICVP16p, label='NICTAL', color='red')
    s3.set_xlabel('Time [min]')
    s3.set_ylabel('Concentration [nM]')
    s3.legend(loc="upper right")
    lg = s3.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)
    
    f4.set_facecolor('white')    
    s4=f4.add_subplot(111)
    s4.plot(TIME, NEPm, linestyle='dashed', label='NEPTUN mRNA', color='blue')
    s4.plot(TIME, NEPp, label='NEPTUN', color='blue')
    s4.plot(TIME, MCTm, linestyle='dashed', label='mCITRIN mRNA', color='green')
    s4.plot(TIME, MCTp, label='mCITRIN', color='green')
    s4.plot(TIME, RFPm, linestyle='dashed', label='RFP mRNA', color='red')
    s4.plot(TIME, RFPp, label='RFP', color='red')
    s4.set_xlabel('Time [min]')
    s4.set_ylabel('Concentration [nM]')
    lg = s4.legend(loc="upper right")
    lg.get_frame().set_alpha(0.5)

    
    if self.gDEBUG:
      path = './out/'+self.__class__.__name__
      paths = [path+'_s0_'+str(NONCE)+'.png', path+'_s1_'+str(NONCE)+'.png', path+'_s2_'+str(NONCE)+'.png', path+'_s3_'+str(NONCE)+'.png', path+'_s4_'+str(NONCE)+'.png']
      f0.savefig(paths[0], format ='png')
      f1.savefig(paths[1], format ='png')
      f2.savefig(paths[2], format ='png')
      f3.savefig(paths[3], format ='png')
      f4.savefig(paths[4], format ='png')
    
    else:
      path = '/var/www/imodel/img/'+self.__class__.__name__
      paths = [path+'_s0_'+str(NONCE)+'.png', path+'_s1_'+str(NONCE)+'.png', path+'_s2_'+str(NONCE)+'.png', path+'_s3_'+str(NONCE)+'.png', path+'_s4_'+str(NONCE)+'.png']
      canvas=FigureCanvas(f0)
      f0.savefig(paths[0], format ='png')
      canvas=FigureCanvas(f1)
      f1.savefig(paths[1], format ='png')
      canvas=FigureCanvas(f2)
      f2.savefig(paths[2], format ='png')
      canvas=FigureCanvas(f3)
      f3.savefig(paths[3], format ='png')
      canvas=FigureCanvas(f4)
      f4.savefig(paths[4], format ='png')
    
    return paths

  def run(self):
    TIME = range(self.gTIME)
    #Gene dosage plasmid [plasmids/cell]
    #Relative number to cell 
    #TODO: to use with measurements transfer function
    #TODO: HEK 293
    #Batard 2001
    #*Assumption: uniform distribution
    #*Tunable
    G1 = self.gG1
    G2 = self.gG2
    G3 = self.gG3
    G4 = self.gG4
    G5 = self.gG5
    G6 = self.gG5
    G7 = self.gG5


    #Experimental TAL Effector rates
    #RLA Measurements
    if self.gDEBUG:
      tal97krab = Effector('TAL97:KRAB', 'dat/tal97krab.dat')
      tal57krab = Effector('TAL57:KRAB', 'dat/tal57krab.dat')
      tal97vp16 = Effector('TAL97:KRAB', 'dat/tal95vp16_4.dat')
      tal57vp16 = Effector('TAL57:KRAB', 'dat/tal95vp16_4.dat')
      nictalkrab = Effector('NICTAL:KRAB', 'dat/tal97krab.dat')
      nictalvp16 = Effector('NICTAL:KRAB', 'dat/tal95vp16_4.dat')
      tal95krab = Effector('TAL57:KRAB', 'dat/tal57krab.dat')
      
      
    else:
      tal97krab = Effector('TAL97:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal97krab.dat')
      tal57krab = Effector('TAL57:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal57krab.dat')
      tal97vp16 = Effector('TAL97:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal95vp16_4.dat')
      tal57vp16 = Effector('TAL57:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal95vp16_4.dat')
      nictalkrab = Effector('NICTAL:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal57krab.dat')
      nictalvp16 = Effector('NICTAL:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal95vp16_4.dat')
      tal95krab = Effector('TAL57:KRAB', '/home/gem/imodel/cellulator/igem-cellulator/Measurements/dat/tal57krab.dat')
      
    #tal97krab.plotTransferFunction()
    #tal57krab.plotTransferFunction()
    #show()
    #quit()
    r97krab = (G1+G3+G4+G5)/(G2)
    r57krab = (G2+G3+G4+G6)/(G1)
    r97vp16 = (G2+G6)/(G2)
    r57vp16 = (G1+G5)/(G1)
    rnickrab = (G1+G2+G5+G6)/(G3)
    rnicvp16 = (G3+G4)/(G3)
    r95krab = (G4)/(G5+G6)
    
    qnic = 1.0/(nictalkrab.getTransferRate2(r97krab)+1e-23)
    q95 = 1.0/(tal95krab.getTransferRate2(r97krab)+1e-23)
    q97 = 1.0/(tal97krab.getTransferRate2(r97krab)+1e-23)
    q57 = 1.0/(tal57krab.getTransferRate2(r57krab)+1e-23)
    #a97 = tal97vp16.getTransferRate2(r97vp16)/1000.0
    #a57 = tal57vp16.getTransferRate2(r57vp16)/1000.0
    a97 = 0.001
    a57 = 0.001
    anic = 0.001
    
    '''
    tal95krab.plotTransferFunction()
    tal97krab.plotTransferFunction()
    tal57krab.plotTransferFunction()
    nictalkrab.plotTransferFunction()
    tal97vp16.plotTransferFunction()
    tal57vp16.plotTransferFunction()
    nictalvp16.plotTransferFunction()
    print q97
    print q57
    print a97
    print a57
    show()
    quit()
    '''



    #Inducer affinity
    #(testing)
    #TODO: look for source
    Kpc = 1.0
    Krm = 1.0

    #Maximal transcription rates (mRNA/min)
    #apply pCMV, pMin
    #TODO: Jan!
    #Apply protein length
    pmin = 1
    pil  = 1
    k1 = 30.0*pmin
    k2 = 30.0*pmin
    k3 = 30.0
    k4 = 30.0*pmin
    k5 = 30.0
    k6 = 30.0
    k7 = 30.0


    #Basal activity of inducible promoter
    #TODO: source, RLA, mammalian
    #TODO: Jan, Miha
    alfapc = 0.085
    alfarm = 0.085
    #alfarm = 0.056

    #Binding affinity 
    #Literature data. PIP repressor.
    #Tigges.
    KPC = 3.0
    KRM = 3.0

    #Cooperativity
    #TODO: source
    npc = 2.0
    nrm = 2.0

    #Translation rates (protein/min)
    #TODO: source
    #TODO: protein length adjust
    t57 = 0.02
    t95 = 0.02
    t97 = 0.02
    tnic = 0.02
    tpc = 0.02
    trm = 0.02
    tnep = 0.02
    tmct = 0.02
    trfp = 0.02


    #Degradation rate (mRNA/min)
    #TODO: source
    kdm = 0.0173

    #Degradation rate (protein/min)
    #TODO: source
    kdp = 0.0058

    
    #Degradation rate (inducer/min)
    kdpc = 2/(exp(self.gpc_half_life))
    kdrm = 2/(exp(self.grm_half_life))
    
    #Concentration scaling factor for copies/cell
    #TODO: source
    fv = 1.85e-3*0.5


    #Temporary switch
    #TODO: sources
    r1 = k1*fv
    r2 = k2*fv
    r3 = k3*fv
    r4 = k4*fv
    r5 = k5*fv
    r6 = k5*fv
    r7 = k5*fv

  

    '''
    OBSERVED SPECIES
    '''
    HACTm = array([0.0 for i in TIME]) #rapamicin binding protein mRNA
    HACTp = array([0.0 for i in TIME]) #rapamicin binding protein
    HACTa = array([0.0 for i in TIME]) #rapamicin binding protein (active)

    PBSm = array([0.0 for i in TIME]) #pristinamicin binding protein mRNA
    PBSp = array([0.0 for i in TIME]) #pristinamicin binding protein
    PBSa = array([0.0 for i in TIME]) #pristinamicin binding protein (active)

    PC = array([0.0 for i in TIME]) #pristinamicin
    RM = array([0.0 for i in TIME]) #rapamicin

    RM[self.grm_start] = 1500
    PC[self.gpc_start] = 1500
            
    T97KRABm = array([0.0 for i in TIME])
    T97KRABp = array([0.0 for i in TIME])
    T95KRABm = array([0.0 for i in TIME])
    T95KRABp = array([0.0 for i in TIME])
    T57KRABm = array([0.0 for i in TIME])
    T57KRABp = array([0.0 for i in TIME])
    TNICKRABm = array([0.0 for i in TIME])
    TNICKRABp = array([0.0 for i in TIME])
    
    T97VP16m = array([0.0 for i in TIME])
    T97VP16p = array([0.0 for i in TIME])
    T57VP16m = array([0.0 for i in TIME])
    T57VP16p = array([0.0 for i in TIME])
    TNICVP16m = array([0.0 for i in TIME])
    TNICVP16p = array([0.0 for i in TIME])
    
    NEPm = array([0.0 for i in TIME])
    NEPp = array([0.0 for i in TIME])
    MCTm = array([0.0 for i in TIME])
    MCTp = array([0.0 for i in TIME])
    RFPm = array([0.0 for i in TIME])
    RFPp = array([0.0 for i in TIME])
    
    
    '''
    SIMULATION
    '''

    
    for t in TIME:
      #Skip first and las step
      if t == 0 or t==self.gTIME-1:
        continue


      #INDUCERS
      #TODO: check with article
      dPBSm = G7 * r7 - kdm * PBSm[t-1]
      dPBSp = tpc * PBSm[t-1] - kdp * PBSp[t-1]
      dPBSa = PBSp[t-1] * (PC[t-1]/(Kpc + PC[t-1]))
      #TODO: check with article
      dHACTm = G7 * r7 - kdm * HACTm[t-1]
      dHACTp = trm * HACTm[t-1] - kdp * HACTp[t-1]
      dHACTa = HACTp[t-1] * (RM[t-1]/(Krm + RM[t-1]))
      r4 = fv * k4 * (alfapc + 120*(1-alfapc)*(PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc))*(HACTa[t-1]**nrm/(KRM + HACTa[t-1]**nrm)))
      r5 = fv * k5 * (alfapc + 120*(1-alfapc)*(PBSa[t-1]**npc/(KPC + PBSa[t-1]**npc)))
      r6 = fv * k6 * (alfarm + 120*(1-alfarm)*(HACTa[t-1]**nrm/(KRM + HACTa[t-1]**nrm)))
      #print 't:'+str(t), r4, r5, r6
      
   
   
      #G1 * r1 - q57 * T57KRABp[t-1] + a57*T57VP16p[t-1]
      #G2 * r2 - q97 * T97KRABp[t-1] + a97*T97VP16p[t-1]
      #G3 * r3 - qnic * TNICKRABp[t-1] + anic*TNICVP16p[t-1]
      #G4 * r4
      #G5 * r5 - q95 * T95KRABp[t-1]
      #G6 * r6 - q95 * T95KRABp[t-1]
      #G7 * r7
      
      
      #TAL97:KRAB
      #r97krab = (G1+G3+G4+G5)/(G2)
      dT97KRABm   = G1 * r1 - q57 * T57KRABp[t-1] + a57*T57VP16p[t-1]
      T97KRABm[t] = T97KRABm[t-1] + dT97KRABm
      T97KRABm[t] = int(T97KRABm[t] > 0) * T97KRABm[t]
      dT97KRABm   = G3 * r3 - qnic * TNICKRABp[t-1] + anic*TNICVP16p[t-1]
      T97KRABm[t] = T97KRABm[t] + dT97KRABm
      T97KRABm[t] = int(T97KRABm[t] > 0) * T97KRABm[t]
      dT97KRABm   = G4 * r4
      T97KRABm[t] = T97KRABm[t] + dT97KRABm
      T97KRABm[t] = int(T97KRABm[t] > 0) * T97KRABm[t]
      dT97KRABm   = G5 * r5 - q95 * T95KRABp[t-1]
      T97KRABm[t] = T97KRABm[t] + dT97KRABm
      T97KRABm[t] = int(T97KRABm[t] > 0) * T97KRABm[t]
      dT97KRABm   = -kdm * T97KRABm[t-1]
      T97KRABm[t] = T97KRABm[t] + dT97KRABm
      T97KRABm[t] = int(T97KRABm[t] > 0) * T97KRABm[t]
      dT97KRABp   = t97 * T97KRABm[t-1] - kdp * T97KRABp[t-1]
      

      
      #TAL57:KRAB
      #r57krab = (G2+G3+G4+G6)/(G1)
      dT57KRABm   = G2 * r2 - q97 * T97KRABp[t-1] + a97*T97VP16p[t-1]
      T57KRABm[t] = T57KRABm[t-1] + dT57KRABm
      T57KRABm[t] = int(T57KRABm[t] > 0) * T57KRABm[t]
      dT57KRABm   = G3 * r3 - qnic * TNICKRABp[t-1] + anic*TNICVP16p[t-1]
      T57KRABm[t] = T57KRABm[t-1] + dT57KRABm
      T57KRABm[t] = int(T57KRABm[t] > 0) * T57KRABm[t]
      dT57KRABm   = G4 * r4
      T57KRABm[t] = T57KRABm[t-1] + dT57KRABm
      T57KRABm[t] = int(T57KRABm[t] > 0) * T57KRABm[t]
      dT57KRABm   = G6 * r6 - q95 * T95KRABp[t-1]
      T57KRABm[t] = T57KRABm[t-1] + dT57KRABm
      T57KRABm[t] = int(T57KRABm[t] > 0) * T57KRABm[t]
      dT57KRABm   = -kdm * T57KRABm[t-1]
      T57KRABm[t] = T57KRABm[t] + dT57KRABm
      T57KRABm[t] = int(T57KRABm[t] > 0) * T57KRABm[t]
      dT57KRABp   = t57 * T57KRABm[t-1] - kdp * T57KRABp[t-1]
      
      
      #TAL97:VP16
      #r97vp16 = (G2+G6)/(G2)
      dT97VP16m = G2 * r2 - q97 * T97KRABp[t-1] + a97*T97VP16p[t-1]
      T97VP16m[t] = T97VP16m[t-1] + dT97VP16m
      T97VP16m[t] = int(T97VP16m[t] > 0) * T97VP16m[t]
      dT97VP16m = G6 * r6 - q95 * T95KRABp[t-1]
      T97VP16m[t] = T97VP16m[t] + dT97VP16m
      T97VP16m[t] = int(T97VP16m[t] > 0) * T97VP16m[t]
      dT97VP16m =  - kdm * T97VP16m[t-1]
      T97VP16m[t] = T97VP16m[t] + dT97VP16m
      T97VP16m[t] = int(T97VP16m[t] > 0) * T97VP16m[t]
      dT97VP16p = t97 * T97VP16m[t-1] - kdp * T97VP16p[t-1]
      
      
      #TAL57:VP16
      #r57vp16 = (G1+G5)/(G1)
      dT57VP16m = G1 * r1 - q57 * T57KRABp[t-1] + a57*T57VP16p[t-1]
      T57VP16m[t] = T57VP16m[t-1] + dT57VP16m
      T57VP16m[t] = int(T57VP16m[t] > 0) * T57VP16m[t]
      dT57VP16m = G5 * r5 - q95 * T95KRABp[t-1]
      T57VP16m[t] = T57VP16m[t-1] + dT57VP16m
      T57VP16m[t] = int(T57VP16m[t] > 0) * T57VP16m[t]
      dT57VP16m   = - kdm * T57VP16m[t-1]
      T57VP16m[t] = T57VP16m[t] + dT57VP16m
      T57VP16m[t] = int(T57VP16m[t] > 0) * T57VP16m[t]
      dT57VP16p = t57 * T57VP16m[t-1] - kdp * T57VP16p[t-1]

      
            
      #NICTAL:KRAB
      #rnickrab = (G1+G2+G5+G6)/(G3)
      dTNICKRABm   = G1 * r1 - q57 * T57KRABp[t-1] + a57*T57VP16p[t-1]
      print dTNICKRABm
      TNICKRABm[t] = TNICKRABm[t-1] + dTNICKRABm
      TNICKRABm[t] = int(TNICKRABm[t] > 0) * TNICKRABm[t]
      dTNICKRABm   = G2 * r2 - q97 * T97KRABp[t-1] + a97*T97VP16p[t-1]
      print dTNICKRABm
      TNICKRABm[t] = TNICKRABm[t] + dTNICKRABm
      TNICKRABm[t] = int(TNICKRABm[t] > 0) * TNICKRABm[t]
      dTNICKRABm   = G5 * r5 - q95 * T95KRABp[t-1]
      print dTNICKRABm
      TNICKRABm[t] = TNICKRABm[t] + dTNICKRABm
      TNICKRABm[t] = int(TNICKRABm[t] > 0) * TNICKRABm[t]
      dTNICKRABm   = G6 * r6 - q95 * T95KRABp[t-1]
      print dTNICKRABm
      TNICKRABm[t] = TNICKRABm[t] + dTNICKRABm
      TNICKRABm[t] = int(TNICKRABm[t] > 0) * TNICKRABm[t]
      dTNICKRABm   = -kdm * TNICKRABm[t-1]
      print dTNICKRABm
      TNICKRABm[t] = TNICKRABm[t] + dTNICKRABm
      TNICKRABm[t] = int(TNICKRABm[t] > 0) * TNICKRABm[t]
      dTNICKRABp   = tnic * TNICKRABm[t-1] - kdp * TNICKRABp[t-1]
      
      
      #NICTAL:VP16
      #rnicvp16 = (G3+G4)/(G3)
      dTNICVP16m   = G3 * r3 - qnic * TNICKRABp[t-1] + anic*TNICVP16p[t-1]
      TNICVP16m[t] = TNICVP16m[t-1] + dTNICVP16m
      TNICVP16m[t] = int(TNICVP16m[t] > 0) * TNICVP16m[t]
      dTNICVP16m   = G4 * r4
      TNICVP16m[t] = TNICVP16m[t] + dTNICVP16m
      TNICVP16m[t] = int(TNICVP16m[t] > 0) * TNICVP16m[t]
      dTNICVP16m   =  - kdm * TNICVP16m[t-1]
      TNICVP16m[t] = TNICVP16m[t] + dTNICVP16m
      TNICVP16m[t] = int(TNICVP16m[t] > 0) * TNICVP16m[t]
      dTNICVP16p   = tnic * TNICVP16m[t-1] - kdp * TNICVP16p[t-1]
      
      
      #TAL95:KRAB
      #r95krab = (G4)/(G5+G6)
      dT95KRABm   = G4 * r4
      T95KRABm[t] = T95KRABm[t-1] + dT95KRABm
      T95KRABm[t] = int(T95KRABm[t] > 0) * T95KRABm[t]
      dT95KRABm   = -kdm * T95KRABm[t-1]
      T95KRABm[t] = T95KRABm[t] + dT95KRABm
      T95KRABm[t] = int(T95KRABm[t] > 0) * T95KRABm[t]
      dT95KRABp   = t95 * T95KRABm[t-1] - kdp * T95KRABp[t-1]
      
      #NEPTUN
      dNEPm = G1 * r1 - q57 * T57KRABp[t-1] - kdm * NEPm[t-1] + a57*T57VP16p[t-1]
      dNEPp = tnep * NEPm[t-1] - kdp * NEPp[t-1]
      
      #MCITRIN
      dMCTm = G2 * r2 - q97 * T97KRABp[t-1] - kdm * MCTm[t-1] + a57*T97VP16p[t-1]
      dMCTp = tmct * MCTm[t-1] - kdp * MCTp[t-1]

      #RFP
      dRFPm = G3 * r3 - qnic * TNICKRABp[t-1] -kdm * RFPm[t-1] + anic*TNICVP16p[t-1] 
      dRFPp = trfp * RFPm[t-1] - kdp * RFPp[t-1]

      #PRISTINAMICYN
      dPC = -PC[t-1] * kdpc
      #RAPAMICIN
      dRM = -RM[t-1] * kdrm
      
      
      
      
      PBSm[t] = PBSm[t-1] + dPBSm
      PBSp[t] = PBSp[t-1] + dPBSp
      HACTm[t] = HACTm[t-1] + dHACTm
      HACTp[t] = HACTp[t-1] + dHACTp   
      PBSa[t] = PBSa[t] + dPBSa
      HACTa[t] = HACTa[t] + dHACTa
      T97KRABp[t] = T97KRABp[t-1] + dT97KRABp
      T97VP16p[t] = T97VP16p[t-1] + dT97VP16p
      T57VP16p[t] = T57VP16p[t-1] + dT57VP16p
      T57KRABp[t] = T57KRABp[t-1] + dT57KRABp
      TNICKRABp[t] = TNICKRABp[t-1] + dTNICKRABp
      TNICVP16p[t] = TNICVP16p[t-1] + dTNICVP16p
      T95KRABp[t] = T95KRABp[t-1] + dT95KRABp
      NEPm[t] = NEPm[t-1] + dNEPm
      NEPp[t] = NEPp[t-1] + dNEPp
      MCTm[t] = MCTm[t-1] + dMCTm
      MCTp[t] = MCTp[t-1] + dMCTp  
      RFPm[t] = RFPm[t-1] + dRFPm
      RFPp[t] = RFPp[t-1] + dRFPp
      
      if not t in self.gpc_start:
        PC[t] = PC[t-1] + dPC
      if not t in self.grm_start:
        RM[t] = RM[t-1] + dRM

      #Preserve above zero
      PBSm[t] = int(PBSm[t] > 0) * PBSm[t]
      PBSp[t] = int(PBSp[t] > 0) * PBSp[t]
      HACTm[t] = int(HACTm[t] > 0) * HACTm[t]
      HACTp[t] = int(HACTp[t] > 0) * HACTp[t]
      PBSa[t] = int(PBSa[t] > 0) * PBSa[t]
      HACTa[t] = int(HACTa[t] > 0) * HACTa[t]
      T97KRABp[t] = int(T97KRABp[t] > 0) * T97KRABp[t]
      T97VP16p[t] = int(T97VP16p[t] > 0) * T97VP16p[t]
      T57KRABp[t] = int(T57KRABp[t] > 0) * T57KRABp[t]
      T57VP16p[t] = int(T57VP16p[t] > 0) * T57VP16p[t]
      TNICKRABp[t] = int(TNICKRABp[t] > 0) * TNICKRABp[t]
      TNICVP16p[t] = int(TNICVP16p[t] > 0) * TNICVP16p[t]
      T95KRABp[t] = int(T95KRABp[t] > 0) * T95KRABp[t]
      NEPm[t] = int(NEPm[t] > 0) * NEPm[t]
      NEPp[t] = int(NEPp[t] > 0) * NEPp[t]
      MCTm[t] = int(MCTm[t] > 0) * MCTm[t]
      MCTp[t] = int(MCTp[t] > 0) * MCTp[t]
      RFPm[t] = int(RFPm[t] > 0) * RFPm[t]
      RFPp[t] = int(RFPp[t] > 0) * RFPp[t]
      PC[t] = int(PC[t] > 0) * PC[t]
      RM[t] = int(RM[t] > 0) * RM[t]

      print T97KRABm[t]
      
    return self.plotResults(TIME, T97KRABm, T97KRABp, T57KRABm, T57KRABp, T97VP16m, T97VP16p, T57VP16m, T57VP16p,  TNICKRABm, TNICKRABp, TNICVP16m, TNICVP16p, T95KRABm, T95KRABp, PBSm, PBSp, PBSa, HACTm, HACTp, HACTa, PC, RM, NEPm, NEPp, MCTm, MCTp, RFPm, RFPp)
    


      

